/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package xenon3d.vector;

import java.io.Serializable;

/**
 * A single precision floating point 4 by 4 matrix.
 * Primarily to support 3D transformations.
 * @author Volker Everts
 * @version 0.1 - 18.09.2011: Created
 */
public class Matrix4f implements Serializable {

    // <editor-fold defaultstate="collapsed" desc=" Static Attributes ">

    /** The serial version UID. */
    private static final long serialVersionUID = 20110918L;

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Private Fields ">

    /** The first matrix element in the first row. */
    float m11 = 1.0f;

    /** The second matrix element in the first row. */
    float m12 = 0.0f;

    /** The third matrix element in the first row. */
    float m13 = 0.0f;

    /** The fourth matrix element in the first row. */
    float m14 = 0.0f;

    /** The first matrix element in the second row. */
    float m21 = 0.0f;

    /** The second matrix element in the second row. */
    float m22 = 1.0f;

    /** The third matrix element in the second row. */
    float m23 = 0.0f;

    /** The fourth matrix element in the second row. */
    float m24 = 0.0f;

    /** The first matrix element in the third row. */
    float m31 = 0.0f;

    /** The second matrix element in the third row. */
    float m32 = 0.0f;

    /** The third matrix element in the third row. */
    float m33 = 1.0f;

    /** The fourth matrix element in the third row. */
    float m34 = 0.0f;

    /** The first matrix element in the fourth row. */
    float m41 = 0.0f;

    /** The second matrix element in the fourth row. */
    float m42 = 0.0f;

    /** The third matrix element in the fourth row. */
    float m43 = 1.0f;

    /** The fourth matrix element in the fourth row. */
    float m44 = 0.0f;

    /** The uniform scale factor. */
    float scale = 1.0f;

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Initialization ">

    /**
     * Constructs and initializes a Matrix4f to the unity matrix.
     */
    public Matrix4f() {

    }

    /**
     * Creates a new Matrix4f with the same values as the specified other
     * matrix
     * @param mat the other matrix
     */
    public Matrix4f(Matrix4f mat) {
        set(mat);
    }

    /**
     * Creates a new Matrix4f from the supplied float array, starting at the
     * specified offset. The number of data elements expected is 17 (1 element
     * for the scale, 16 for the matrix).<p>
     * At the first position in the data array is the uniform scale factor of
     * the matrix, followed by 16 data elements representing the matrix itself.<p>
     * The matrix data elements are expected in column-major order.<p>
     * NOTE: it is strongly advised to check whether the matrix data really
     * represent a 3D transformation matrix with uniform scale factor. This can
     * be checked by calling isTransformation().
     * @param fv the float array containing the matrix data in column first order
     * @param offset the offset at which the data ist stored in the float array
     */
    public Matrix4f(float[] fv, int offset) {
        set(fv, offset);
    }

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Public Properties (General) ">

    /**
     * Sets the value of this Matrix4f to the value of another matrix.
     * @param mat the other matrix
     */
    public final void set(Matrix4f mat) {
        m11 = mat.m11;
        m12 = mat.m12;
        m13 = mat.m13;
        m14 = mat.m14;
        m21 = mat.m21;
        m22 = mat.m22;
        m23 = mat.m23;
        m24 = mat.m24;
        m31 = mat.m31;
        m32 = mat.m32;
        m33 = mat.m33;
        m34 = mat.m34;
        m41 = mat.m41;
        m42 = mat.m42;
        m43 = mat.m43;
        m44 = mat.m44;
        scale = mat.scale;
    }

    /**
     * Sets the values in this Matrix4f from the supplied float array, starting
     * at the specified offset.<p>
     * Only matrix rows 1 to 3 are expected, because row 4 of a transformation
     * matrix is always (0,0,0,1) and need not to be saved or restored. The
     * number of data elements expected is thus 13 (one element for the scale,
     * followed by 12 for matrix rows 1 to 3).<p>
     * The matrix data elements are expected in column-major order.<p>
     * NOTE: it is strongly advised to check whether the matrix data really
     * represents a 3D transformation matrix with uniform scale factor. This can
     * be checked by calling isTransformation().
     * @param fv the float array containing the matrix data in column first order
     * @param offset the offset at which the data ist stored in the float array
     */
    public final void set(float[] fv, int offset) {
        scale = fv[offset];
        m11 = fv[offset + 1];
        m21 = fv[offset + 2];
        m31 = fv[offset + 3];
        m41 = 0;
        m12 = fv[offset + 4];
        m22 = fv[offset + 5];
        m32 = fv[offset + 6];
        m42 = 0;
        m13 = fv[offset + 7];
        m23 = fv[offset + 8];
        m33 = fv[offset + 9];
        m43 = 0;
        m14 = fv[offset + 10];
        m24 = fv[offset + 11];
        m34 = fv[offset + 12];
        m44 = 1;
    }

    /**
     * Copies the values of this Matrix4f into the supplied float array,
     * starting at the specified offset.<p>
     * Only matrix rows 1 to 3 are written, because row 4 of a transformation
     * matrix is always (0,0,0,1) and need not to be saved or restored. The
     * number of data elements written is thus 13 (one element for the scale,
     * followed by 12 for matrix rows 1 to 3).<p>
     * The matrix data elements are written in column-major order.
     * @param fv the float array that will contain the values of the matrix
     * @param offset the offset at which the data ist to be stored in the float array
     */
    public final void get(float[] fv, int offset) {
        fv[offset + 0] = scale;
        fv[offset + 1] = m11;
        fv[offset + 2] = m21;
        fv[offset + 3] = m31;
        fv[offset + 4] = m12;
        fv[offset + 5] = m22;
        fv[offset + 6] = m32;
        fv[offset + 7] = m13;
        fv[offset + 8] = m23;
        fv[offset + 9] = m33;
        fv[offset + 10] = m14;
        fv[offset + 11] = m24;
        fv[offset + 12] = m34;
    }

    /**
     * Computes the determinate of this matrix.
     * @return the determinate of the matrix
     */
    public final float det() {
        float
            det =  m11 * (m22 * m33 * m44 + m23 * m34 * m42 + m24 * m32 * m43 - m24 * m33 * m42 - m22 * m34 * m43 - m23 * m32 * m44);
            det -= m12 * (m21 * m33 * m44 + m23 * m34 * m41 + m24 * m31 * m43 - m24 * m33 * m41 - m21 * m34 * m43 - m23 * m31 * m44);
            det += m13 * (m21 * m32 * m44 + m22 * m34 * m41 + m24 * m31 * m42 - m24 * m32 * m41 - m21 * m34 * m42 - m22 * m31 * m44);
            det -= m14 * (m21 * m32 * m43 + m22 * m33 * m41 + m23 * m31 * m42 - m23 * m32 * m41 - m21 * m33 * m42 - m22 * m31 * m43);
        return det;
    }

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Public Properties (Identity and Scale) ">

    /**
     * Sets this Matrix4f to identity.
     */
    public final void setIdentity() {
        m11 = 1.0f;
        m12 = 0.0f;
        m13 = 0.0f;
        m14 = 0.0f;
        m21 = 0.0f;
        m22 = 1.0f;
        m23 = 0.0f;
        m24 = 0.0f;
        m31 = 0.0f;
        m32 = 0.0f;
        m33 = 1.0f;
        m34 = 0.0f;
        m41 = 0.0f;
        m42 = 0.0f;
        m43 = 0.0f;
        m44 = 0.0f;
        scale = 1.0f;
    }

    /**
     * Sets the value of this matrix to a scale matrix with the passed scale.
     * Any preexisting rotation and translation will be reset to 0.
     * @param scale the scale factor for the matrix
     */
    public final void set(float scale) {
        setIdentity();
        m11 = scale;
        m22 = scale;
	m33 = scale;
        this.scale = scale;
    }

    /**
     * Sets the scale component of the current matrix by factoring out the
     * current scale and multiplying by the new scale. Any rotational and
     * translational components are preserved.
     * @param scale the new scale
     */
    public final void setScale(float scale) {
        if (scale < FMath.EPS) throw new IllegalArgumentException();
        float s = scale / this.scale;
        m11 *= s;
        m12 *= s;
        m13 *= s;
        m21 *= s;
        m22 *= s;
        m23 *= s;
        m31 *= s;
        m32 *= s;
        m33 *= s;
        this.scale = scale;
    }

    /**
     * Returns the uniform scale factor.
     * @return the uniform scale factor of this matrix
     */
    public final float getScale() {
        return scale;
    }

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Public Properties (Position) ">

    /**
     * Sets the value of this matrix to a position transform matrix with the
     * passed position. Any preexisting rotation will be reset to 0, the scale
     * will be reset to 1.
     * @param p the position transform point
     */
    public final void set(Point3f p) {
        setIdentity();
        m14 = p.x;
        m24 = p.y;
        m34 = p.z;
    }

    /**
     * Sets the translational component of this matrix. Any rotational and scale
     * components are preserved.
     * @param x the new x coordinate
     * @param y the new x coordinate
     * @param z the new x coordinate
     */
    public final void setPosition(float x, float y, float z) {
        m14 = x;
        m24 = y;
        m34 = z;
    }

    /**
     * Sets the translational component of this matrix. Any rotational and scale
     * components are preserved.
     * @param pos the new position
     */
    public final void setPosition(Point3f pos) {
        m14 = pos.x;
        m24 = pos.y;
        m34 = pos.z;
    }

    /**
     * Returns the translational component of this matrix.
     * @param result the point that will receive the translational component
     */
    public final void getPosition(Point3f result) {
        result.x = m14;
        result.y = m24;
        result.z = m34;
    }

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Public Properties (Rotation - Absolute) ">

    /**
     * Sets the value of this matrix to a counter clockwise rotation about the
     * x axis. All of the non-rotational components are set as if this were an
     * identity matrix.
     * @param angle the angle in degrees to rotate about the X axis
     */
    public final void rotX(float angle) {
        float sin = FMath.sin(angle);
        float cos = FMath.cos(angle);
        setIdentity();
        m22 = cos;
        m23 = -sin;
        m32 = sin;
        m33 = cos;
    }

    /**
     * Sets the value of this matrix to a counter clockwise rotation about the
     * y axis. All non-rotational components are reset as if this were an
     * identity matrix.
     * @param angle the angle to rotate about the Y axis
     */
    public final void rotY(float angle) {
        float sin = FMath.sin(angle);
        float cos = FMath.cos(angle);
        setIdentity();
        m11 = cos;
        m13 = sin;
        m31 = -sin;
        m33 = cos;
    }

    /**
     * Sets the value of this matrix to a counter clockwise rotation about the
     * z axis. All non-rotational components are reset as if this were an
     * identity matrix.
     * @param angle the angle to rotate about the Z axis
     */
    public final void rotZ(float angle) {
        float sin = FMath.sin(angle);
        float cos = FMath.cos(angle);
        setIdentity();
        m11 = cos;
        m12 = -sin;
        m21 = sin;
        m22 = cos;
    }

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Public Properties (Rotation - EulerAngles) ">

    /**
     * Sets the value of this matrix to the rotation matrix converted from the
     * Euler angles provided; Any preexisting scale as well as the translational
     * component is preserved.
     * The euler rotation angles are applied first about the X, then Y then Z
     * axis. These rotations are applied using a static frame of reference.
     * In other words, the orientation of the Y rotation axis is not affected
     * by the X rotation and the orientation of the Z rotation axis is not
     * affected by the X or Y rotation.
     * @param e the EulerAngle3f
     */
    public final void setRotation(EulerAngle3f e) {
        float cr = FMath.cos(e.x);
        float sr = FMath.sin(e.x);
        float cp = FMath.cos(e.y);
        float sp = FMath.sin(e.y);
        float cy = FMath.cos(e.z);
        float sy = FMath.sin(e.z);
        float srsp = sr * sp;
        float crsp = cr * sp;
        m11 = cp * cy;
        m12 = cp * sy;
        m13 = -sp;
        m21 = srsp * cy - cr * sy;
        m22 = srsp * sy + cr * cy;
        m23 = sr * cp;
        m31 = crsp * cy + sr * sy;
        m32 = crsp * sy - sr * cy;
        m33 = cr * cp;
        if (scale != 1.0f) {
            float s = scale;
            scale = 1.0f;
            setScale(s);
        }
    }

    /**
     * Returns an euler rotation tuple that is equivalent to that set by
     * <code>setRotation(Tuple3f)</code>. Note that it does not necessarily
     * return the *same* Euler angles, but the rotation will be equivalent, i.e.
     * will have the same result when used to rotate a vector or node.
     * @param result the euler angles that are equivalent to the rotation
     * represented by this matrix.
     */
    public final void getRotation(EulerAngle3f result) {
        result.set(this);
    }

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Public Properties (Rotation - AxisAngle) ">

    /**
     * Sets the value of this matrix to the matrix conversion of the (single
     * precision) axis and angle arguments. Any preexisting scale is as well as
     * the translational component is preserved.
     * @param angle the angle in degrees
     * @param axis the axis vector
     */
    public final void setRotation(Vector3f axis, float angle) {
        float mag = axis.length();
        if (mag < FMath.EPS) setIdentity();
        else {
            float ax = axis.x / mag;
            float ay = axis.y / mag;
            float az = axis.z / mag;
            float sin = FMath.sin(angle);
            float cos = FMath.cos(angle);
            float t = 1.0f - cos;
            float xz = ax * az;
            float xy = ax * ay;
            float yz = ay * az;
            m11 = t * ax * ax + cos;
            m12 = t * xy - sin * az;
            m13 = t * xz + sin * ay;
            m21 = t * xy + sin * az;
            m22 = t * ay * ay + cos;
            m23 = t * yz - sin * ax;
            m31 = t * xz - sin * ay;
            m32 = t * yz + sin * ax;
            m33 = t * az * az + cos;
        }
        if (scale != 1.0f) {
            float s = scale;
            scale = 1.0f;
            setScale(s);
        }
    }

    /**
     * Gets the axis of rotation and an angle that are equivalent to the axis
     * and angle pair set by <code>setRotation(float, Vector3f)</code>. Note
     * that this method does not necessarily return the *same* angle and axis,
     * but the rotation will be equivalent, i.e. will have the same result when
     * used to rotate a vector or node.
     * @param axis the axis of rotation
     * @return the angle of rotation in degrees
     */
    public final float getRotation(Vector3f axis) {
        float x = m32 - m23;
        float y = m13 - m31;
        float z = m21 - m12;
        float angle = 0.0f;
        float mag = FMath.sqrt(x * x + y * y + z * z);
        if (mag > FMath.EPS) {
            float sin = 0.5f * mag;
            float cos = 0.5f * (m11 + m22 + m33 - 1.0f);
            axis.x = x / mag;
            axis.y = y / mag;
            axis.z = z / mag;
            angle = FMath.atan2(sin, cos);
        }
        else axis.set(0.0f, 1.0f, 0.0f);
        return angle;
    }

    /**
     * Sets the value of this matrix to the matrix conversion of the (single
     * precision) axis-angle. Any preexisting scale as well as the translational
     * component is preserved.
     * @param a the axis-angle object
     */
    public final void setRotation(AxisAngle4f a) {
        Vector3f axis = new Vector3f();
        float angle = a.get(axis);
        setRotation(axis, angle);
    }

    /**
     * Gets the axis-angle which is equivalent to the axis-angle set by
     * <code>setRotation(AxisAngle4f)</code>. Note that this method does not
     * necessarily return the *same* axis-angle, but the rotation will be
     * equivalent, i.e. will have the same result when used to rotate a vector
     * or node.
     * @param result the axis-angle
     */
    public final void getRotation(AxisAngle4f result) {
        result.set(this);
    }

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Public Properties (Rotation - Quaternion ">

    /**
     * Sets the value of this matrix to the matrix conversion of the (single
     * precision) quaternion argument. Any preexisting scale as well as the
     * translational component is preserved.
     * @param q the quaternion to be converted
     */
    public final void setRotation(Quat4f q) {
        float x = q.x;
        float y = q.y;
        float z = q.z;
        float w = q.w;
        m11 = 1.0f - 2.0f * y * y - 2.0f * z * z;
        m12 = 2.0f * (x * y - w * z);
        m13 = 2.0f * (x * z + w * y);
        m21 = 2.0f * (x * y + w * z);
        m22 = 1.0f - 2.0f * x * x - 2.0f * z * z;
        m23 = 2.0f * (y * z - w * x);
        m31 = 2.0f * (x * z - w * y);
        m32 = 2.0f * (y * z + w * x);
        m33 = 1.0f - 2.0f * x * x - 2.0f * y * y;
        if (scale != 1.0f) {
            float s = scale;
            scale = 1.0f;
            setScale(s);
        }
    }

    /**
     * Returns the normalized rotational component of this matrix; the values
     * are placed into the Quat4f parameter.
     * @param result quaternion into which the rotation component is placed
     */
    public final void getRotation(Quat4f result) {
        result.set(this);
    }

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Public Properties (Rotation - Matrix) ">

    /**
     * Sets the normalized rotational component of this matrix. Any preexisting
     * scale as well as the translational component is preserved.
     * @param mat matrix which contains the rotational component
     */
    public final void setRotation(Matrix3f mat) {
        float s = scale;
        setRotationScale(mat);
        if (scale != s) setScale(s);
    }

    /**
     * Returns the normalized rotational component of this matrix; the values
     * are placed into the Matrix3d parameter.
     * @param result matrix into which the rotational component is placed
     */
    public final void getRotation(Matrix3f result) {
        getRotationScale(result);
        result.normalize();
    }

    /**
     * Replaces the upper 3x3 matrix values of this matrix with the values in
     * the specified matrix. The translational component is preserved.
     * @param mat the matrix that will be the new upper 3x3
     */
    public final void setRotationScale(Matrix3f mat) {
        m11 = mat.m11;
        m12 = mat.m12;
        m13 = mat.m13;
        m21 = mat.m21;
        m22 = mat.m22;
        m23 = mat.m23;
        m31 = mat.m31;
        m32 = mat.m32;
        m33 = mat.m33;
        scale = mat.scale;
    }

    /**
     * Returns the upper 3x3 values of this matrix and places them into the
     * result matrix
     * @param result the matrix that will hold the values
     */
    public final void getRotationScale(Matrix3f result) {
        result.m11 = m11;
        result.m12 = m12;
        result.m13 = m13;
        result.m21 = m21;
        result.m22 = m22;
        result.m23 = m23;
        result.m31 = m31;
        result.m32 = m32;
        result.m33 = m33;
        result.scale = scale;
    }

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Public Properties (Transformable - AxisAngle) ">

    /**
     * Sets the value of this matrix from a rotation expressed by an axis-angle
     * and a position. The scale is set to 1.
     * @param rot the rotation expressed as an axis-angle quaternion
     * @param pos the position expressed as a point
     */
    public final void setTransform(AxisAngle4f rot, Point3f pos) {
        setRotation(rot);
        setPosition(pos);
    }

    /**
     * Sets the value of this matrix from a rotation expressed by an axis-angle
     * a position, and a uniform scale factor.
     * @param rot the rotation expressed as an axis-angle
     * @param pos the position expressed as a point
     * @param s the scale value
     */
    public final void setTransform(AxisAngle4f rot, Point3f pos, float s) {
        setRotation(rot);
        setPosition(pos);
        setScale(s);
    }

    /**
     * Returns the normalized rotational component of this matrix, the
     * translation, and the scale. None of the matrix values are modified.
     * @param rot the normalized result axis-angle representing the rotation
     * @param pos the translation result component
     * @return the uniform scale factor of this matrix
     */
    public final float getTransform(AxisAngle4f rot, Point3f pos) {
        getRotation(rot);
        getPosition(pos);
        return scale;
    }

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Public Properties (Transformable - Quaternion) ">

    /**
     * Sets the value of this matrix from the rotation expressed by the
     * quaternion q and the translation point p. The scale is set to 1.
     * @param rot the rotation expressed as a quaternion
     * @param pos the position expressed as a point
     */
    public final void setTransform(Quat4f rot, Point3f pos) {
        setRotation(rot);
        setPosition(pos);
    }

    /**
     * Sets the value of this matrix from the rotation expressed by the
     * quaternion q, the translation t, and the scale s.
     * @param rot the rotation expressed as a quaternion
     * @param pos the position expressed as a point
     * @param s the scale value
     */
    public final void setTransform(Quat4f rot, Point3f pos, float s) {
        setRotation(rot);
        setPosition(pos);
        setScale(s);
    }

    /**
     * Returns the normalized rotational component of this matrix, the
     * translation, and the scale. None of the matrix values are modified.
     * @param rot the normalized result quaternion representing the rotation
     * @param pos the translation result component
     * @return the scale component of this transform
     */
    public final float getTransform(Quat4f rot, Point3f pos) {
        getRotation(rot);
        getPosition(pos);
        return scale;
    }

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Public Properties (Transformable - Matrix) ">

    /**
     * Sets the rotational component of this matrix (perhaps containing a scale
     * factor), and the specified translation.
     * @param rot the matrix representing the rotation and scale
     * @param pos the translation component expressed as a point
     */
    public final void setTransform(Matrix3f rot, Point3f pos) {
        setRotationScale(rot);
        setPosition(pos);
    }

    /**
     * Sets the normalized rotational component of this matrix, the
     * translation, and the scale.
     * @param rot the normalized result matrix representing the rotation
     * @param pos the translation component expressed as a point
     * @param scale the scale component of this transform
     */
    public final void setTransform(Matrix3f rot, Point3f pos, float scale) {
        setRotation(rot);
        setPosition(pos);
        setScale(scale);
    }

    /**
     * Returns the scaled rotational component of this matrix, the translation,
     * and the uniform scale factor. None of the matrix values are modified.
     * @param rot the result matrix representing the rotation
     * @param pos the translation result component
     * @return the scale component of this matrix
     */
    public final float getTransform(Matrix3f rot, Point3f pos) {
        getRotation(rot);
        getPosition(pos);
        return scale;
    }

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Public Methods (Transformation / Normalization ">

    /**
     * Performs singular value decomposition normalization of the rotational
     * component of this matrix. Resets the scale to 1, but preserves any
     * existing rotation and position.
     */
    public final void normalize() {
        if (scale != 1.0f) {
            m11 /= scale;
            m12 /= scale;
            m13 /= scale;
            m21 /= scale;
            m22 /= scale;
            m23 /= scale;
            m31 /= scale;
            m32 /= scale;
            m33 /= scale;
            scale = 1.0f;
        }
        m41 = 0.0f;
        m42 = 0.0f;
        m43 = 0.0f;
        m44 = 1.0f;
    }

    /**
     * Perform cross product normalization of this matrix. Resets the scale to 1,
     * but preserves any existing rotation and position.
     */
    public final void normalizeCP() {
        float mag = 1.0f / FMath.sqrt(m11 * m11 + m21 * m21 + m31 * m31);
        m11 *= mag;
        m21 *= mag;
        m31 *= mag;
        mag = 1.0f / FMath.sqrt(m12 * m12 + m22 * m22 + m32 * m32);
        m12 *= mag;
        m22 *= mag;
        m32 *= mag;
        m13 = m21 * m32 - m22 * m31;
        m23 = m12 * m31 - m11 * m32;
        m33 = m11 * m22 - m12 * m21;
        scale = 1.0f;
        m41 = 0.0f;
        m42 = 0.0f;
        m43 = 0.0f;
        m44 = 1.0f;
    }

    /**
     * Transforms the specified point p by multiplying this matrix by p and
     * places the result back into p. (p = this * p).
     * @param p the point to be multiplied by this matrix and then replaced
     */
    public final void transform(Point3f p) {
        transform(p, p);
    }

    /**
     * Transforms the specified point p by multiplying this matrix by p and and
     * places the result into another point, the result point. (result = this * p).
     * @param p  the point to be multiplied by this matrix
     * @param result  the point into which the transformed product is placed
     */
    public final void transform(Point3f p, Point3f result) {
        float x = p.x;
        float y = p.y;
        float z = p.z;
        result.x = m11 * x + m12 * y + m13 * z + m14;
        result.y = m21 * x + m22 * y + m23 * z + m24;
        result.z = m31 * x + m32 * y + m33 * z + m34;
    }

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Public Methods (Matrix Math) ">

    /**
     * Sets the value of this matrix to its transpose.
     */
    public final void transpose() {
        float temp;
        temp = m21;
        m21 = m12;
        m12 = temp;
        temp = m31;
        m31 = m13;
        m13 = temp;
        temp = m41;
        m41 = m14;
        m14 = temp;
        temp = m32;
        m32 = m23;
        m23 = temp;
        temp = m42;
        m42 = m24;
        m24 = temp;
        temp = m43;
        m43 = m34;
        m34 = temp;
    }

    /**
     * Inverts this matrix in place.
     */
    public final void invert() {
        transpose();
    }

    /**
     * Sets the value of this matrix to the result of multiplying itself with
     * another matrix.
     * @param mat the other matrix
     */
    public final void mul(Matrix4f mat) {
        mul(this, mat);
    }

    /**
     * Sets the value of this matrix to the result of multiplying the two
     * argument matrices together.
     * @param m1 the first matrix
     * @param m2 the second matrix
     */
    public final void mul(Matrix4f m1, Matrix4f m2) {
        float a01 = m1.m11 * m2.m11 + m1.m12 * m2.m21 + m1.m13 * m2.m31 + m1.m14 * m2.m41;
        float a02 = m1.m11 * m2.m12 + m1.m12 * m2.m22 + m1.m13 * m2.m32 + m1.m14 * m2.m42;
        float a03 = m1.m11 * m2.m13 + m1.m12 * m2.m23 + m1.m13 * m2.m33 + m1.m14 * m2.m43;
        float a04 = m1.m11 * m2.m14 + m1.m12 * m2.m24 + m1.m13 * m2.m34 + m1.m14 * m2.m44;
        float a05 = m1.m21 * m2.m11 + m1.m22 * m2.m21 + m1.m23 * m2.m31 + m1.m24 * m2.m41;
        float a06 = m1.m21 * m2.m12 + m1.m22 * m2.m22 + m1.m23 * m2.m32 + m1.m24 * m2.m42;
        float a07 = m1.m21 * m2.m13 + m1.m22 * m2.m23 + m1.m23 * m2.m33 + m1.m24 * m2.m43;
        float a08 = m1.m21 * m2.m14 + m1.m22 * m2.m24 + m1.m23 * m2.m34 + m1.m24 * m2.m44;
        float a09 = m1.m31 * m2.m11 + m1.m32 * m2.m21 + m1.m33 * m2.m31 + m1.m34 * m2.m41;
        float a10 = m1.m31 * m2.m12 + m1.m32 * m2.m22 + m1.m33 * m2.m32 + m1.m34 * m2.m42;
        float a11 = m1.m31 * m2.m13 + m1.m32 * m2.m23 + m1.m33 * m2.m33 + m1.m34 * m2.m43;
        float a12 = m1.m31 * m2.m14 + m1.m32 * m2.m24 + m1.m33 * m2.m34 + m1.m34 * m2.m44;
        float a13 = m1.m41 * m2.m11 + m1.m42 * m2.m21 + m1.m43 * m2.m31 + m1.m44 * m2.m41;
        float a14 = m1.m41 * m2.m12 + m1.m42 * m2.m22 + m1.m43 * m2.m32 + m1.m44 * m2.m42;
        float a15 = m1.m41 * m2.m13 + m1.m42 * m2.m23 + m1.m43 * m2.m33 + m1.m44 * m2.m43;
        float a16 = m1.m41 * m2.m14 + m1.m42 * m2.m24 + m1.m43 * m2.m34 + m1.m44 * m2.m44;
        m11 = a01;
        m12 = a02;
        m13 = a03;
        m14 = a04;
        m21 = a05;
        m22 = a06;
        m23 = a07;
        m24 = a08;
        m31 = a09;
        m32 = a10;
        m33 = a11;
        m34 = a12;
        m41 = a13;
        m42 = a14;
        m43 = a15;
        m44 = a16;
    }

    /**
     * Returns true if the L-infinite distance between this matrix and another
     * matrix m1 is less than or equal to the specified positive epsilon
     * parameter. The L-infinite distance is equal to MAX[i=0,1,2 ; j=0,1,2 ;
     * abs(this.m(i,j) - mat.m(i,j)].
     * @param mat the matrix to be compared to this matrix
     * @param eps the epsilon parameter
     * @return true, if this matrix is epsilon equal to the other matrix
     */
    public boolean equals(Matrix4f mat, float eps) {
        return
            FMath.isEqual(m11, mat.m11, eps) &&
            FMath.isEqual(m12, mat.m12, eps) &&
            FMath.isEqual(m13, mat.m13, eps) &&
            FMath.isEqual(m14, mat.m14, eps) &&
            FMath.isEqual(m21, mat.m21, eps) &&
            FMath.isEqual(m22, mat.m22, eps) &&
            FMath.isEqual(m23, mat.m23, eps) &&
            FMath.isEqual(m24, mat.m24, eps) &&
            FMath.isEqual(m31, mat.m31, eps) &&
            FMath.isEqual(m32, mat.m32, eps) &&
            FMath.isEqual(m33, mat.m33, eps) &&
            FMath.isEqual(m34, mat.m34, eps) &&
            FMath.isEqual(m41, mat.m41, eps) &&
            FMath.isEqual(m42, mat.m42, eps) &&
            FMath.isEqual(m43, mat.m43, eps) &&
            FMath.isEqual(m44, mat.m44, eps);
    }

    // </editor-fold>

    // <editor-fold defaultstate="collapsed" desc=" Implementation Object ">

    @Override
    public int hashCode() {
        int bits = 1;
        bits = 31 * bits + Float.floatToIntBits(m11);
        bits = 31 * bits + Float.floatToIntBits(m12);
        bits = 31 * bits + Float.floatToIntBits(m13);
        bits = 31 * bits + Float.floatToIntBits(m14);
        bits = 31 * bits + Float.floatToIntBits(m21);
        bits = 31 * bits + Float.floatToIntBits(m22);
        bits = 31 * bits + Float.floatToIntBits(m23);
        bits = 31 * bits + Float.floatToIntBits(m24);
        bits = 31 * bits + Float.floatToIntBits(m31);
        bits = 31 * bits + Float.floatToIntBits(m32);
        bits = 31 * bits + Float.floatToIntBits(m33);
        bits = 31 * bits + Float.floatToIntBits(m34);
        bits = 31 * bits + Float.floatToIntBits(m41);
        bits = 31 * bits + Float.floatToIntBits(m42);
        bits = 31 * bits + Float.floatToIntBits(m43);
        bits = 31 * bits + Float.floatToIntBits(m44);
        return bits;
    }

    @Override
    public boolean equals(Object obj) {
        if (this == obj) return true;
        if (obj == null) return false;
        if (obj instanceof Matrix4f) {
            Matrix4f mat = (Matrix4f) obj;
            return
                m11 == mat.m11 && m12 == mat.m12 && m13 == mat.m13 && m14 == mat.m14 &&
                m21 == mat.m21 && m22 == mat.m22 && m23 == mat.m23 && m24 == mat.m24 &&
                m31 == mat.m31 && m32 == mat.m32 && m33 == mat.m33 && m34 == mat.m34 &&
                m41 == mat.m41 && m42 == mat.m42 && m43 == mat.m43 && m44 == mat.m44;
        }
        return false;
    }

    /**
     * Returns a string that contains the values of this Matrix3f.
     * @return the String representation
     */
    @Override
    public String toString() {
        return
            "[" + m11 + ", " + m12 + ", " + m13 + ", " + m14 + "]\n" +
            "[" + m21 + ", " + m22 + ", " + m23 + ", " + m24 + "]\n" +
            "[" + m31 + ", " + m32 + ", " + m33 + ", " + m34 + "]\n" +
            "[" + m41 + ", " + m42 + ", " + m43 + ", " + m44 + "]";
    }

    // </editor-fold>

} // end class Matrix4f